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9 Abstract. In a companion paper, equations for partially molten media 

10 were derived using two-scale homogenization theory. This approach begins 

11 with a grain scale description and then coarsens it through multiple scale ex- 

12 pansions into a macroscopic model. One advantage of homogenization is that 

13 effective material properties, such as permeability and the shear and bulk 

« viscosity of the two-phase medium, are characterized by cell problems, bound- 
is ary value problems posed on a representative microstructural cell. The so- 
le lutions of these problems can be averaged to obtain macroscopic parameters 
I? that are consistent with a given microstructure. This is particularly impor- 

18 tant for estimating the "compaction length" which depends on the product 

19 of permeability and bulk viscosity and is the intrinsic length scale for vis- 

20 cously deformable two-phase flow. 

21 In this paper, we numerically solve ensembles of cell problems for several 

22 geometries. We begin with simple intersecting tubes as this is a one param- 

23 eter family of problems with well known results for permeability. Using this 

24 data, we estimate relationships between the porosity and all of the effective 

25 parameters with curve fitting. For this problem, permeability scales as <p n , 

26 n ~ 2 — 3, as expected and the bulk viscosity scales as <j)~ m , m ~ 1, which 

27 has been speculated, but never shown directly for deformable porous media. 

28 The second set of cell problems add spherical inclusions at the tube inter- 

29 sections. These show that the permeability is controlled by the smallest pore 

30 throats and not by the total porosity, as expected. The bulk viscosity remains 
si inversely proportional to the porosity, and we conjecture that this quantity 
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32 is insensitive to the specific microstructure. The computational machinery 

33 developed can be applied to more general geometries, such as texturally equi- 

34 librated pore shapes. However, we suspect that the qualitative behavior of 

35 our simplified models persists in these more realistic structures. In partic- 

36 ular, our hybrid numerical-analytical model predicts that for purely mechan- 

37 ical coupling at the microscale, all homogenized models will have a compaction 

38 length that vanishes as porosity goes to zero. This has implications for com- 

39 putational models, and it suggests that these models might not resist com- 
« plete compaction. 
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1. Introduction 

41 Partially molten regions in the Earth's asthenosphere (e.g. beneath mid-ocean ridges or 

42 subduction zones) are usually modeled as a viscously deformable permeable media. Such 

43 models are typically composed of macroscopic equations for the conservation of mass, 

44 momentum and energy of each phase. In our companion paper, Simpson et al. [2008a], 

45 we derived several systems of governing equations for partially molten systems using 

46 homogenization. Briefly, we began with a grain scale description of two interpenetrating 

47 fluids, each satisfying the Stokes equations, coupled by their common interface. Several 
4a different macroscopic models could then be coarsened from this microscopic description, 
49 depending on our assumptions on the velocities, viscosities, and grain scale geometry. 

so An important feature of this approach is that macroscopic properties such as perme- 

51 ability, shear viscosity, and bulk viscosity naturally appear in the macroscopic equations, 

52 even if they are not defined at the grain scale. In particular, permeability and bulk vis- 

53 cosity are properties of the two-phase aggregate, not the volume averages of small scale 

54 variations. In contrast, previous work on the magma problem, including McKenzie [1984]; 

55 Bercovici and Ricard [2003, 2005]; Bercovici [2007]; Hier-Majumder et al. [2006]; Ricard 
se [2007], began with models much larger than the grain scale. There, the viscosities, per- 
s? meability, and other closures were assumed and justified from other results. In contrast, 
sa homogenization derives these properties self-consistently. 

59 While homogenization techniques appropriately inserts the constitutive relationships 

60 into the macroscopic equations, connecting them to the microstructure requires the solu- 

61 tion of specific "cell problems." For the physical system derived in Simpson et al. [2008a], 
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there are actually a series of ten Stokes problems for both fluid and solid posed on the 
micro-scale (which can be reduced to four for micro-structures with sufficient symmetry). 

In this paper, we numerically explore the cell problems to extract parameterizations 
of the various effective parameters in terms of porosity, a simple measurement of the 
microstructure. A feature of this work is to derive constitutive relationships for the bulk 
viscosity as a function of porosity, which is essential for describing compactible permeable 
media. In particular, we show that for a range of simple pore structure the effective bulk 
viscosity is related to the porosity as. 

Ccff. OC <f X 

Our results and some additional theory suggest that this scaling is insensitive to the 
specific pore geometry. 

Since permeability and bulk-viscosity can both be consistently related to the same mi- 
crostructure, these calculations also allow us to study the "compaction length" [McKenzie, 
1984]. At sufficiently low melt concentrations, it is approximately 

<5comp. OC \/k e s.(eS. 

which depends on the product of the derived permeability k c s, and bulk viscosity Ceff., 
both of which depend on the porosity. The compaction length is the intrinsic length scale 
in magma dynamics. This work suggests that under solely mechanical deformation, 

lim5 comp .(0) = 0, 

<f>— >o 

implying that no mechanical mechanism prevents a region from compacting to zero poros- 
ity. This result also places strong resolution constraints on computational models of 
magma migration which may require a regularization for small porosities. 
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eg An outline of this work is as follows. In Section 2, we review several models of partially 

70 molten rock and highlight the constitutive relations. In Section 3, we demonstrate the 

?! process of assuming a cell geometry to extract computational closures for the macroscopic 

72 system. We revisit these closures in Section 4 for more general geometries to assess their 

73 robustness. Finally, in Section 5, we combine our numerics with the equations and examine 

74 the implications. 

2. Review of Equations and Constitutive Relations 
2.1. Macroscopic Equations 

In Simpson et al. [2008a], we showcased three models for momentum conservation in a 

partially molten medium. They were distinguished by the assumed scalings for the relative 

velocities and viscosities between the fluid and solid phases, along with the connectivity 

of the pore network. One of them, dubbed Biphasic-I, was given by the equations: 



75 Notation for his model may be found in Table 1. In particular, /c eff and Ceff. are the 

76 emergent permeability and bulk viscosity. i] cS . is an auxiliary, tensorial, shear viscosity 

77 capturing grain scale anisotropy. All are related to the aforementioned cell problems. 

78 We highlight this case because (la - lc) is nearly identical to the models of McKenzie, 

79 Bercovici, Ricard, and others in the absence of melting and surface physics. 

2.2. Constitutive Relations 
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so The constitutive relations for the permeability and viscosity are fundamental to the 

si dynamics of these models. Indeed, they are the source of much nonlinearity and it is 

82 useful to review some proposed closures. Notation for this appears in Table 2. 

ss At low porosity, it is common to relate permeability, k, to porosity by a power law, 

84 k oc 0™. Estimates of n vary, n ~ 2 — 5, Scheidegger [1974]; Bear [1988]; Dullien [1992]; 

8 5 Turcotte and Schubert [2002]; McKenzie [1984]; Doyen [1988]; Cheadle [1989]; Martys 

86 et al. [1994]; Faul et al. [1994]; Faul [1997, 2000]; Koponen et al. [1997]; Wark and Watson 

87 [1998]; Wark et al. [2003]; Cheadle et al. [2004]. For partially molten rocks, the exponent 

88 is better constrained by both analysis and experiment to n ~ 2 — 3. 

For the matrix shear viscosity, Hirth and Kohlstedt [1995a, b]; Kohlstedt et al. [2000]; 
Kelemen et al. [1997]; Kohlstedt [2007] experimentally observed a melt weakening effect, 
which they fit to the curve: 

H 9+ f = fis exp (-0/0*) , 0* = O(10~ 2 ) (2) 

89 fj, s+ f is the shear viscosity of the solid matrix in the presence of melt; /i s is the shear 

90 viscosity a melt-free matrix. In Bercovici et al. [2001]; Bercovici and Ricard [2003] and 

91 related works, the viscosity is weighted by (1 — 0), which is also present in (la). Reiterat- 

92 ing, fi s+ f is a fitting of experimental data. Regardless, the shear viscosity is taken to be 

93 isotropic, and porosity weakening in other models. 

94 Lastly, the bulk viscosity, ( s , is often taken as ( s oc _m , m ~ — 1, though m is usually 

95 either zero or one. m = has often been used because the variation of bulk viscosity 

96 with porosity is poorly constrained by observations. Scott and Stevenson [1984] and oth- 

97 ers invoked the bore hole studies of ice by Nye [1953] to justify m — 1. In that work, 

98 Nye considered the dynamics of individual spherical and cylindrical voids in an infinite 
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99 medium. Taylor [1954]; Prud'homme and Bird [1978] computed m — 1, in the limit of 

wo small porosity, using models of incompressible fluids mixed with gas bubbles. Schmeling 

101 [2000] also employed ( s oc citing studies for the analogous question of the effective 

102 bulk modulus of a fluid filled poroelastic medium. Those results rely on the self-consistent 

103 approximation methodology, treated in Torquato [2002]. Models in Bercovici et al. [2001]; 

104 Bercovici and Ricard [2003] and related works possess a term that functions as a 

105 bulk viscosity, though it has a very different origin. McKenzie [1984] suggested using a 
we metallurgical model of spheres due to Arzt et al. [1983] which considered a range of mech- 
10? anisms for densification of powders under hot isostatic pressing. For viscous compaction 

108 mechanisms, the Arzt result recover a bulk viscosity proportional to _1 for Newtonian 

109 viscosities. However, for very small porosities, they suggest that surface diffusion effects 
no become important and imply that ( s oc log(0 _1 ). Finally, Connolly and Podladchikov 
m [1998] invoke an assymetric bulk viscosity that is weaker during expansion than during 
n2 compaction, which is motivated by the deformation of brittle crustal materials. However, 
us it is unclear if this rheology is relevant to high-temperature/high-pressure creeping ma- 
n4 terials as are expected in the mantle. In general, it is expected that the bulk-viscosity 
us should become unbounded as the porosity reduces to zero as the system then becomes 
no incompressible. 

2.3. Cell Problems 

n7 In homogenization, a medium with fine scale features is modeled by introducing two or 

us more spatial scales. As discussed in Simpson et al. [2008a], the direct approach makes 

n9 multiple scale expansions of the dependent variables, letting them depend on both the 

120 coarse and fine length scales. 
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The two characteristic lengths in our model are L, the macroscopic scale, and £, the 
grain scale. Their ratio, 

i 

e= L 

121 is the parameter in which we perform the multiple scale expansions. In the process of 

122 performing the expansions and matching powers of e, we are mathematically constrained 

123 to solve a collection of auxiliary cell problems posed on the fine scale. Notation for this 

124 setup is provided in Table 3. 

125 It is analytically advantageous to approximate the mixture as a periodic medium; such 

126 a configuration appears in Figure 1. Q, the macroscopic region containing both matrix 
12? and melt is periodically tiled with scaled copies of the cell, scaled to unity in Figure 2. 

128 The cell problems are posed within Y s , the matrix portion of the cell, and Yf, the melt 

129 portion for the cell. Y s and Yf meet on interface 7. 

Generically, the cell problems take the form: 

V y ■ (-pi + 2e„(v)) = f in y> or Y s (3a) 
V y -v = g in Y f or Y s (3b) 
(—pi + 2e y (v)) • n = r • n or u = U on 7 (3c) 

130 which are generally, compressible Stokes flow problems for the microscopic solid and melt 

131 velocities and pressures (v,p). Here, V y is the gradient, V y - is the divergence, and e y 

132 the strain-rate operator defined on the fine length scale y. A more complete description 

133 is given in Simpson et al. [2008a]. f, g, r, and U are prescribed forcing functions on the 

134 relevant portion of Y, either Y s or Yf. The solution, (v,p), is periodic on the portion of 

135 the boundary not intersecting the interface 7. The cell problems may be interpreted as 
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136 the fine scale response to an applied stress on a unit cell of either the melt or the matrix. 

137 The material parameters - /c e fr., Ccff.) and r] ef[ - are then defined as the cell average of an 
us appropriate manipulation of these pressures and velocities defined at the pore scale. In 
us particular, an average of the solid pressure in one problem determines the bulk viscosity 

140 and certain averages of melt velocities in another problem determine the permeability. 

141 Again, we emphasize that the material parameters are not volume averages of the same 
«2 parameters at the fine scale. 

3. Effective Parameters: Intersecting Tube Geometry 

«3 To connect the effective material properties - i] cS , k cS , and Ccfr. ~ with the porosity, 
«4 we must solve the cell problems and extract a parameterization. Unfortunately, the so- 
ws lution of the Stokes equations in a generic three dimensional domain lacks an analytic 
«6 representation. Thus we compute the solutions to the cell problems numerically and fit 
147 the results to appropriate parameterized constitutive models. Notation for this section is 
us summarized in Table 4. 

149 As a first example of numerically closing the constitutive relations in a homogenization 

150 based model, we study the cell domain of triply intersecting cylinders, pictured in Figure 

151 3. The fluid occupies the cylinders while the solid is the complementary portion of the 

152 cube. While this geometry is an oversimplification of real pore-geometries, it serves as 

153 proof of concept of the unified, self-consistent homogenization algorithm. Additionally, it 

154 serves as a numerical check since the permeability for this problem is expected to scale as 

155 2 . In Section 4 we examine a generalization of this geometry. 
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156 In what follows, we emphasize that the information given by the parameterizations 

157 should be read with a certain skepticism. The order of magnitude and signs of the coeffi- 

158 cients and exponents are of greater utility than the particular numbers. 

We remark here that for the intersecting tube geometry, the tube radius scaled to the 
cell-size, b, can be explicitly related to the porosity: 

= 3nb 2 - 8V2b 3 (4) 

159 We solve our problems using finite elements on unstructured meshes using the open 
wo source libraries from the FEniCS and PETSc projects [e.g. Dupont et al, 2003; Kirby and 

161 Logg, 2007; Logg, 2007; Balay et al, 2004, 2001], with meshes generated using CUBIT 

162 [Sandia Corporation, 2008]. Details of this method and numerical benchmarks are given 

163 in Appendix B. 

3.1. Effective Permeability 

The first cell problem we treat is for permeability. The equations are: 

-V y q l + Vjk 1 = -e, in Y f (5a) 
V y ■ k* = in Y f (5b) 
k* = on 7 (5c) 

where k* is a three-dimensional velocity and q l is a scalar pressure, for each i = 1,2,3. 
These are equations (64a - 64c) in Simpson et al. [2008a] and describe the motion of the 
melt through the porous matrix. In general, the permeability is the second order tensor: 

(K) f = \j Yf Vdy J Yf k 2 dy J Yf k s dy 
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The symmetry properties of the domain simplify this to: 

(K)f = (k})7, k* is the first component of vector k 1 

Thus, it is sufficient to compute the case % = 1. Then /c e ff., the permeability of the matrix, 
in (lb) is 



kcE. = (k\) f \ (6) 

Darcy's Law and permeability have been studied by many techniques, including ho- 
mogenization; we refer the reader to the references in Section 2.2. We study it here to 
understand how the permeability behaves in concert with the other constitutive relations 
as the microstructure varies. This also serves as a benchmark problem for our software; 
see Appendix B. 

As noted, porosity and permeability are often related by a power law, k oc <p n , with 
n ~ 2—5. To motivate such a relation, we turn to a toy model, as presented in Turcotte and 
Schubert [2002]. The melt is assumed to be in Poiseuille flow through triply intersecting 
cylinders. Additionally, the cylinders have small radii; it is a low porosity model. The 
permeability of such a system is 

Ktoy-i = ^ ~ O.OO44£ 2 2 . (7) 

Other simple models are developed in Scheidegger [1974]; Bear [1988]; Dullien [1992]. 

We now fit our computed permeabilities, to porosity by such a relation. For the 
tube domains, the least squares fit is 

(k\} f = exp(-4.42 ± .1O5)0 2 - 2O± - 0391 . (8) 

This curve and the data appear in Figure 4. The fit matches expectations of an O(10~ 3 — 
10~ 2 ) prefactor and an exponent ~ 2 — 3. The error in (8) is the associated 95% confidence 
DRAFT September 4, 2009, 4:18pm DRAFT 
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172 interval. We report these intervals in all regressions, though they rely on the specious 

173 assumption that error in our synthetic data is normally distributed. 

3.2. Effective Bulk Viscosity 

The effective bulk viscosity is related to the solution (£, () of: 

V,-(-C/ + 2e,(0) =0 inr; (9a) 
V w -£ = 1 mY s (9b) 
(-C/ + 2e,(O)-n = on 7 (9c) 

where £ is a three-dimensional velocity and ( is a scalar pressure. These are equations 
(Bfa - Bfc) in Simpson et al. [2008a] and are associated with the compaction of the 
matrix. The effective bulk viscosity is then defined as 



2 



(10) 



174 The dependence of the effective bulk viscosity of partially molten rock as a function of 

175 porosity is the most poorly constrained of the material properties. This is partly due to 

176 the difficulties in constructing an experiment that will measure it as a function of porosity 

177 independently of the shear viscosity McKenzie [1984]; Kelemen et al. [1997]; Stevenson 

178 and Scott [1991]. As mentioned in Section 2.2, a bulk viscosity oc _1 has often appeared 

179 in the literature. 

Two toy models for the bulk viscosity of an incompressible fluid seeded with compress- 
ible gas bubbles were formulated by Taylor [1954]; Prud'homme and Bird [1978]. They 



DRAFT 



September 4, 2009, 4:18pm 



DRAFT 



X - 14 SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

relate the bulk viscosity to the porosity as: 

Cs = lj Taylor (11a) 

( s = ^ (1 - <p) Prud'homme and Bird (lib) 

wo Taylor's expression, (11a), relied on a single inclusion model for a gas bubble in an infinite 

lai medium, (lib) is derived by considering a sphere of fluid with a gas filled spherical cavity, 

182 and seeking the bulk viscosity of a compressible fluid that will give rise to the same 

183 radial stress for specified boundary motion. The 1 — is due to Prud'homme and Bird 

184 restricting their model to a finite volume. This factor also appears in the proposed bulk 
las viscosity of Schmeling [2000]. These expressions have motivated the use of ( s oc in 
we macroscopic models of partial melts, although they actually arise from a different problem 
is? of a compressible inclusion in an incompressible fluid, rather than the divergent flow of two 
las interconnected incompressible fluids. Surprisingly, when we consider a simple toy problem 
189 to approximate the the full cell calculation, we find that the expressions are identical. 

Appendix A provides details of this toy problem, which is related to equations (9a - 
9c), and gives the solution 

< c >«=§K i+ C) (i - 0) ' 

from which we get 

a = £<w). 

This motivates trying to numerically fit (Q 8 to (1 — 0) p /0 9 , expecting p and q to be close 
to unity. Indeed, the data, plotted in Figure 5, fits the curve 

( C ) s = eX p(-0.131 ± O.OO514)0- LO2±aooi32 (l - 0)O.884±o.oo869 (12) 

wo which is quite similar to the scaling of the toy model. 
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3.3. Supplementary Anisotropic Viscosity 

We now examine the cell problem related to the supplementary viscosity r] e e_, a fourth 
order tensor. The equations are: 

V y ■ (-7T /m J + 2e y (x lm )) =0 in Y s (13a) 
V y -x lm = inY s (13b) 
(-TT^dij + 2e yAj {x lm )) n 3 = ~\ Wjm + S im 6ji) rij on 7 (13c) 

where x lm is a three-dimensional velocity vector and n lm is a scalar pressure, for each 
pair (l,m), I = 1,2,3 and m = 1,2,3. These are equations (B2a - B2c) from Simpson 
et al. [2008a] and are tied to tensorial surface stresses applied on the matrix. Using the 
solutions (x lrn , ir lm ), 

Vll = (e y (x lm ))s (14) 



wi Because of symmetry, we need only consider two problems: (l,m) = (1,1) and (l,m) = 

192 (1,2), corresponding to normal stress and shear stress in each direction. 

Although we have no toy problem as motivation, P (1 — <f)) q proved to be satisfactory. 
First, we study the problem (l,m) = (1, 1), a uniaxial stress problem. For the tubes, we 
fit 

-(eyff)), = exp(-1.72 ± .O4O5)0- 964± - olO4 (l - 0)^.0685 (15) 

193 This vanishes as <fi — > and as <fi — > 1 and it is nearly linear at small porosity. The curves 

194 and the data are plotted in Figure 6. 

There is also the simple shear stress problem, (l,m) = (1,2). For this, we fit 

_(e 12 (^( 12 ))) s = exp(-1.04 ± .OISS)^' ^- 0048 ^! - 0)i-"±.O3i8 (16) 

195 Plots for this are given in Figure 7. 
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we We could now employ the simple computational closures (8), (12), (15), and (16), in 
197 (la - lc), to simulate and study the macroscopic problem. 

4. Generalization of the Cell Domains 

Ms Regrettably, Earth materials are not as trivial as intersecting cylinders. Even an ide- 

199 alized olivine grain is a tetrakaidekahedron, pictured in Figures 8. As depicted, some 

200 fraction of the melt lies along the triple junctions and some is at the quadruple junc- 

201 tions. Other examples of idealized, texturally equilibrated arrangements appear in von 

202 Bargen and Waff [1986] and Cheadle [1989]; Cheadle et al. [2004]. These methods are also 

203 amenable to studying random media. One could compute relations based on ensembles 

204 of randomly generated cell domains, solving all relevant cell problems on the ensemble. 

Motivated by Figure 8, we explore a simple generalization of the tube geometry by 
adding a sphere of independent radius at the intersection, as in Figure 9. This retains 
the symmetry of the previous model, but adds a second parameter, allowing multiple 
geometries for the same porosity. The sphere captures some aspect of the pocket at the 
quadruple junctions. The sphere radius, a, and the tube radius, b, are related to the 
porosity by the equation: 



205 We shall refer to it as the sphere+tube geometry. 

206 We now repeat the computations of Section 3 on the sphere+tube geometry to assess 

207 the sensitivity of the effective parameters to cell geometry. We also perform computations 

208 on domains where the fluid occupies an isolated sphere at the center of a cube. Though 

209 this is disconnected, it provides useful information. 
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4.1. Permeability Revisited 

10 Parameterizing the porosity-permeability relation for this generalization is quite chal- 

11 lenging. In contrast to the tube geometry, the data points of the sphere+tube geometry, 

12 plotted in Figure 4, do not collapse onto a curve. There is some positive correlation be- 

13 tween permeability and porosity, and for a given porosity, the permeability of the equiv- 
« alent tube geometry is an upper bound. 

15 To better understand the trend, we examine the computed flow fields in Figures 10 

16 and 11. These plot the velocity magnitude on two fluid domains with the same tube 

17 size, but different sphere sizes. Most of the flow is within the tube. While there is some 
is detrainment as it enters the sphere, the flow in the tube appears insensitive to the size of 
is the sphere. 

This motivates fitting permeability against tube radius. Indeed, an alternative to (7), 

is 

^toy-II = ^ (18) 

The tube diameter 5, is equivalent to 2b, b the the tube radius in the tube and sphere+tube 
geometries. Both data sets appear in Figure 12. This is a significant improvement over 
Figure 4. The least square fits are: 

(k\) f = exp(-0.592 ± .0354)6 4 - 10± - 0156 , for tube geometry (19) 
(k[) f = exp(-0.628 ± .198)& 3 93± 0684 , for sphere+tube geometry (20) 

20 These estimates with (18); taking 5 = 2b and scaling out £, this relationship is fc t0 y-ii = 

21 .125& 4 . The data is still positively correlated with sphere radius, altering it by as as much 

22 as an order of magnitude. The deviations are greatest when both b <C 1 and b <^ a. 



DRAFT 



September 4, 2009, 4:18pm 



DRAFT 



X - 18 SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

223 That the permeability is more strongly correlated with the tube radius than the overall 

224 geometry is not surprising. Koponen et al. [1997] discuss the notion of effective porosity, 

225 the portion of the void space where there is significant flow. Denoting our effective porosity 

226 cff ., we seek a relation k e g. oc 4>™ s 

Given the flow fields in Figures 10 and 11 and the success with the tube radius fittings, 
we posit that the effective porosity is the portion of the porosity within the tubes. A two- 
dimensional analog appears in Figure 13. Using (8), we define c g. for the sphere+tube 
domains: 

eff = 3tt6 2 - 8V2b 3 (21) 

227 Zhu and Hirth [2003] made a similar approximation; from von Bargen and Waff [1986] 

228 they construed that the permeability was controlled by the minimal cross-sectional area 

229 of the pore network. A similar argument is made by Cheadle [1989]. In our domains, the 

230 minimal cross-sectional area is 7rb 2 . 

We fit 

(k\) f = exp(-4.44 ± .144)</4 06± - 0374 for sphere+tube. (22) 

231 This appears in Figure 14. Again, deviation is highest for very large spheres with very 

232 thin tubes. Unfortunately, e ff. does not satisfy a conservation law, making it a less than 

233 ideal macroscopic quantity to track, though it does satisfy the bound (f> e s, < 0. 

There is still as much as an order of magnitude deviation at low porosity from relation 
(22). The unresolved part of the permeability for the 4> e s. fit is increasing in the sphere 
radius. This motivates trying to fit against both cff and another parameter. It is sufficient 
to fit permeability to c g. and 0, resulting in 

(k}) / = exp(-4.20 ± .O681)0 L ff 88± - O22 V 351± - 0300 for sphere+tube. (23) 
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This is plotted in Figure 15. Some deviation persists at low porosity, but it is less than an 
order of magnitude. Both the sphere+tube data points and the tube data points collapse 
onto this curve. (23) is also consistent with (8), the fit of porosity agaist permeability for 
the tube geometry. Taking eff . = for the tubes, (23) becomes 

(k\) f = exp(-4.2O)0 2 - 23 (24) 

This is similar to the most general permeability relationships, formulated in Scheidegger 
[1974]; Bear [1988]: 

permeability = £ 2 /i(pore shape)/ 2 (0) (25) 

234 By including both <p and c g. in (23), we capture some aspect of the pore shape. 

4.2. Bulk and Supplementary Viscosities Revisited 

In contrast to the permeability problem, the effective bulk viscosity and supplementary 
anisotropic viscosity are quite robust to the domain distortion. As before, we fit (() s to 
(1 — <f)) p /(j) g , expecting p and q to be close to unity. For the sphere+tube and the sphere 
geometry, the least squares fits are: 

(C)s = exp(0.301 ± O.OlO2)0- LOO±aooi74 (l - 0)°- 718±0 - 0337 , for sphere geometry (26a) 

(C)s = exp(0.124 ± O.O975)0- a985±ao252 (l - 0) 1 - O9±O - 186 ! for sphere+tube geometry 

(26b) 

235 The data and these fits are plotted appear in Figure 5. The spherical geometry appears 

236 to be an upper bound on the bulk viscosity for a given porosity. We also remark that the 

237 prefactors vary by less than an order of magnitude amongst the different domains. This 

238 is a strong endorsement of an effective bulk viscosity oc not only for small porosity, 

239 but also for moderate porosities > 10%. 



DRAFT 



September 4, 2009, 4:18pm 



DRAFT 



X - 20 SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

The supplementary anisotropic viscosity terms are also robust under this geometric 
perturbation. For the problem (l,m) = (1, 1), the two new geometries fit: 

-{ei,i(x U ))s = exp(-1.68 ± .O588)0- 98O± - O1O1 (1 - 0) 3 - 56± - 195 , for sphere geometry 

(27a) 

-(ei i(x n )) s = exp(-1.94 ± .139)0- 912± - O359 (1 - 0) L25± - 265 , for sphere+tube geometry 

(27b) 

2« The curves and the data are also plotted in Figure 6. For < 10%, the spread amongst 

241 the three geometries is less than an order of magnitude. 

Similar results are found in the (/, m) = (1, 2) problem. The two new additional domains 
are fit with: 

_(e 12 (^( 12 ))) s = exp(-1.67 ± .O222)0 LO2± - OOO38O (1 - 0)°-4oo±.o7 37) for gphere geometry 

(28a) 

-(ei 2 (x (12) )) s = exp(-1.32 ± .OOSSS)^ 1 - 03 ^ 00228 ^ - ^o-sni.ies^ for sphere+tube geometry 

(28b) 

242 The sphere+tube data is bounded between the sphere data and the tube data; the spread 

243 is less than an order of magnitude. 

5. Discussion and Open Problems 

244 We have successfully parameterized the macroscopic parameters on ensemble of domains 

245 for the simple pore geometries. These can now be consistently used with the macroscopic 

246 model given by equations (la - lc). We now review and discuss our computations, both 

247 independently of and together with the macroscopic equations. 

5.1. Sensitivity to Geometry 

248 As demonstrated by our computations in Sections 3 and 4, the effective parameters 

249 demonstrate a variety of sensitivities to the geometry. Permeability, for our cell geome- 

250 tries, can be bounded by porosity, fc e ff. J$ 0™, where n ~ 2, as seen in Figure 4. But in 
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251 general, it cannot be expressed as a function of a single parameter, such as porosity. How- 

252 ever, if the shapes are constrained by textural equilibration, as in Cheadle [1989], there 

253 will be a single variable parameterization for each dihedral angle. Solving the cell prob- 

254 lems on these shapes, and assessing their sensitivity to the dihedral angle, is an important 

255 open problem. 

256 In contrast, the bulk viscosity appears to be rather insensitive to the geometry, scaling 

257 as Ccfi. oc — 0). We see this from the variation, or lack thereof, between the data 

258 and fits for the isolated spheres and the triply intersecting cylinders, in Figure 5. Though 

259 these are entirely different geometric structures, there is less than an order of magnitude of 

260 variation in the computed Cefr.- These too merit computation on the texturally equilibrated 

261 shapes. Randomly generated geometries may also be of interest. 

262 Both components of the supplementary anisotropic shear viscosity appear to be insen- 

263 sitive to geometry, with less than an order of magnitude amongst the three geometries. 

264 However, because the problems are driven by surface stresses, it may be that a more 

265 anisotropic shape could alter these scalings. Studying them on randomly generated shapes 

266 may provide insight on the role of grain scale anisotropy. 

5.2. Bulk Viscosity 

267 Perhaps our most significant result is the self-consistent bulk viscosity, arising from 

268 a purely mechanical model of partially molten rock. A spatially varying bulk viscosity 

269 is quite important. Indeed, significant differences in dynamics were noted between the 

270 solutions of the McKenzie [1984] model and the model in Ricard et al. [2001]. The authors 

271 point to the use of a constant bulk viscosity in the McKenzie model as the source of the 

272 discrepancy. 
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273 Prior to Ricard et al. [2001], Schmeling [2000] remarked that the _1 dependence has 

274 an important impact on the compaction length. For <p = 0(1%), the bulk viscosity is two 

275 orders of magnitude greater than the shear viscosity. While many studies took ( s and fi s 

276 to be the same order, this higher bulk viscosity leads to a compaction length an order of 

277 magnitude greater. Ricard [2007] made a similar observation on the impact of variable 

278 bulk viscosity on the compaction length. 

279 Schmeling also commented that this variable bulk viscosity could induce melt focusing 

280 towards the axis in his plume simulations. This is an additional nonlinearity that may 

281 be important to geophysical problems. Many studies relying on the McKenzie model 

282 employed a constant bulk viscosity, including Richter and McKenzie [1984]; Spiegelman 

283 and McKenzie [1987]; Spiegelman [1993a, b, c]; Aharonov et al. [1997]; Spiegelman [1996]; 

284 Kelemen et al. [1997]; Spiegelman et al. [2001]; Katz et al. [2004]; Spiegelman et al. [2007]. 

285 Spiegelman and Kelemen [2003]; Spiegelman [2003] used a bulk viscosity, with ( s oc (p~ m 

286 with m > n, n the exponent in the permeability relationship k oc <p n to prevent the system 

287 from compacting to zero between the reactive channels. It would be interesting to revisit 

288 these problems with a 0" 1 bulk viscosity. 

5.3. Compaction Length 

We now combine (la - lc), our leading order equations derived by multiple scale ex- 
pansions, with our computed constitutive relations. For <fi <C 1 our numerical estimates, 



DRAFT 



September 4, 2009, 4:18pm 



DRAFT 



SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS X - 23 

(23), (26a-26b), (15-27b), (27a), and (16-28a) are approximately: 

Ceflf. ~ /i s Co0 _1 (l - 0) 
VeS. ~ /V?O0(1 - 0) 

k is a 0(1CT 3 — 1CT 2 ) constant, Co is an 0(1) constant, and r] is a O(10 _1 ) constant 
fourth order tensor. Under these assumptions, the equations for the Biphasic-I model, 
(la - lc), simplify: 

O = pg-VP + V[/i s Co0^V-V s ] 



+ V- 



2(l-0)/i s e(V s )--(l-0)/i s V-V s J 



+ V- [2// s 0(l - <j>)r^e xM {Y 8 )] 

(v/-v^-MW! (V p_ g /) 

Pi 

V • [0V / + (1 - 0)V S ] = 



(29) 

(30) 
(31) 



If we were to use these equations and numerically derived constitutive relations to solve 
a boundary value problem, the compaction length would again appear as an important 
length scale, 



[Ceflf. + |a*s(1 - 0) + 2 |?7eff.|] fceff. 



'comp. 



'^(1 - 0XCO0- 1 + | + 2 \ Vo \ <P)k <Pli<p^ 



If < 1, then Co0 _1 > 4/3 + 2 1 770 1 and 1-0^1, 



325^.95 
efF. 



(32) 



Since e ff. < 0, we have an upper bound on the compaction length, 



4omp. < h/(okoHs/Hf<fr S 
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Hence, 

lim <5 comp . = 

</>->o 

289 We believe (33), which constrains the compaction length by the porosity raised to a small 

290 positive power, is relatively insensitive to the geometric configuration. This follows from 

291 our alleged robustness of our effective bulk viscosity, Ccfr. oc and the broad agreement 

292 m the porosity-permeability relationship, k oc <p n with n > 2. 

293 A compaction length that vanishes with porosity has interesting consequences. For 

294 example, this compaction length scaling does not rule out the possibility that a partially 

295 molten rock could expel all fluid by mechanical means. It also does not permit the 

296 infiltration of fluid into a dry region. Understanding how any of these systems of equations 

297 transition between a partially molten region and a dry region is an outstanding question. 

298 We also note that though our scaling relationship does not forbid compaction, it does 

299 not imply it either. There may be a dynamic response that prevents the matrix from 

300 mechanically compacting to zero. Such effects were mathematically proven to exist in a 

301 one-dimensional simplification of the model, without melting or freezing, in Simpson et al. 

302 [2007]; Simpson and Weinstein [2008]; Simpson et al. [2008c]. 

303 If, instead, we had concluded <5 comp . oc <p q , with q < 0, then the compaction length 

304 would become unbounded as the melt vanished. Hence the region of deformation in the 

305 matrix needed to segregate additional fluid would also become infinite, precluding further 

306 segregation solely by mechanical processes. 

Because we have parameterized the permeability with e fr., we can explicitly see the re- 
sponse of the compaction length as the melt network becomes disconnected. e ff. measures 
the volume fraction where melt flows. As the channels close up and the melt becomes 
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trapped and e g. — > 0. Taking this limit in (32), 

lim 5 comp . = 

The compaction length can vanish, even if the melt fraction remains bounded away from 
zero. A similar conclusion could be drawn for the compaction length of McKenzie, 



so? Letting k = K (j) n , if we interpret the loss of connectivity as k — > 0, 5ms4 vanishes with 

308 nonzero porosity. 

5.4. Other Physics 

309 There are several results that were not realized by our model, and these merit discus- 

310 sion. We did not recover the ( s oc log(0 _1 ) result of Arzt et al. [1983]. This arises in 

311 the limit of sufficiently low porosity that the primary transport mechanism is by grain 

312 boundary diffusion. Since we did not include any surface physics in our fine scale model, 

313 we should not have expected to upscale their effects. If good grain scale descriptions of 
3« these processes could be formulated, it might be possible to coarsen them via homog- 

315 enization, perhaps recovering this macroscopic relation. However, as log(0" 1 ) becomes 

316 unbounded more slowly than 0~ 1 , this result does not change the basic argument about 
si? the compaction length vanishing with zero porosity. 

Another relationship not captured by either our work is the experimental fit for matrix 
shear viscosity in the presence of melt from Hirth and Kohlstedt [1995a, b]; Kelemen et al. 




(34) 



[1997]; Kohlstedt et al. [2000]; Kohlstedt [2007], 



fx s+f oc exp (-0/0*) 
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318 Our model possesses a porosity weakening mechanism; all of the viscosity terms are oc 

319 1 — 0. However, the anisotropic part, r] eS , is not sign definite and is small compared 

320 to the isotropic component. Furthermore, there does not appear to be an exponential 

321 relation. Hirth and Kohlstedt [1995a] hypothesized that the presence of melt enhances 

322 grain boundary diffusion, providing a fast path for deformation through the melt. As 

323 with the log(0 _1 ) bulk viscosity, this is a surface physics phenomenon not captured by 
32-1 our Stokes models. 

5.5. Open Problems 

325 As discussed, two important open problems are the computation of the cell problems 

326 on more realistic and general cell domains and the inclusion of surface physics into the 

327 model. The former problem is rather straightforward, requiring good computational tools 

328 for generating the domains and solving the Stokes equations on them. The latter problem 

329 is more challenging, requiring fine scale equations for these processes. 

330 The models could also be augmented by giving the matrix a nonlinear rheology, leading 

331 to nonlinear cell problems. Though these could be solved and studied numerically, this 

332 is much more difficult as the different forcing components no longer decouple. Rather 

333 than being able to split the matrix cell problems into a bulk viscosity problem, and two 

33 4 surface stress problems, they would have to be done simultaneously. However, the derived 

335 parameterizations for the effective viscosities would be important to magma migration. 

336 A nonlinear matrix rheology is expected at large strain rates and was needed to compu- 

337 tationally model physical experiments for shear bands in Katz et al. [2006]. 

338 Another important physical rheology is viscoelasticity. Connolly and Podladchikov 

339 [1998]; Vasilyev et al. [1998] extended the earlier, purely viscous, models to include elastic 
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effects. This is an important regime since it permits both short and long time scales, as 
is found at the asthenosphere-lithosphere boundary. 



Appendix A: Spherical Model 



Here we develop a toy model for equations (9a 



9c), 



V, • (-(I + 2e y (0) 







inF s 



i 



inn 



(- C 7 + 2e,(0)- 



n 







on 7 



whose solution yields the effective bulk viscosity, 

2 

Consider a fluid domain, Yf, occupying a small isolated sphere at the center of the unit 
cube; Y s is the complementary region. Smoothing out the exterior boundary of Y s deforms 
it into a sphere. We solve the dilation stress problem on this domain. To avoid confusion, 
let 



342 Our toy problem is posed on Y s sp ere . 

Although periodicity is no longer a meaningful boundary condition, it can be shown 
that the normal velocity on the periodic part of dY s vanishes. We set the normal velocity 
to zero on the exterior boundary of y s s P here . On the interior shell, the stress free condition 



sphere 



{y eR 3 | a< |y| < 1} 



(Al) 



s 




(A2) 
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remains. The equations are: 

V • (-(I + 2e(0) = in Y7 phere (A3a) 

V • f = 1 in F/ phere (A3b) 

(-(I + 2e(0) • n = at r = a (A3c) 

f • n = at r = 1 (A3d) 

a < 1 is the radius of the interior sphere. Decomposing the velocity into incompressible, 
v mc ', and compressible, VII, components, the compressible part solves: 

V 2 n = 1 (A4) 

Let the boundary conditions on the potential be: 

n| r=a = 0, Vn| r= i = (A5) 

Since the problem is spherically symmetric, its solution is 

Ii=\{r 2 -a 2 )+ l -{r- 1 -a- 1 ) (A6) 

The incompressible velocity must be divergence free. Again, by spherical symmetry, 

V • v inc - = ^d r (r 2 4 nc -) = v™- = C/r 2 

To satisfy the boundary condition at r = 1, v mc - • n = t>™ c - = 0. Therefore, C = and 
v™' = 0. The pressure then solves d r ( = 0, so it is constant, 

C(r) = ((a) for r £ (a, 1) 

Applying boundary condition (A3c), 

1 / 4 

C = 2e rr (0|,=a = 2d 2 U\ r=a = - 2 + 



3 V" ' a 3 
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In this geometry, = a 



3 , and the cell averaged pressure is 




(1-0) 



(A7) 



Therefore, 



30 



(1-0) 



(A8) 



343 (A8) uses definition (10), but the £'s in each solve different problems. (A8) is specific to 

344 the geometry of a sphere with cavity, while (10) is for a generic geometry occupying some 

345 fraction of the unit cube. 

346 (A7) is plotted along with data for the numerical solutions of the cell problem posed 

347 on Y s cube in Fi gure 16. There is good agreement between (A7) and these computations 
343 for porosity > 10% suggesting our deformation y s cubc =>. y^sphere wag reasonaD i e . We 

349 reiterate that though our result is quite similar to that of (11a) and (lib), the origin of 

350 the underlying Stokes problem is quite different. 

Appendix B: Computational Methods and Results 
Bl. Cell Problem Boundary Conditions 

Though we could solve the cell problems as stated, with the specified boundary condi- 
tions on interface 7 and periodic on the rest of the domain, we use the symmetry of our 
model problems to reduce the computational cost by a factor of eight. The symmetry 
properties of the solutions allow us to formulate the appropriate boundary conditions on 
the portion of the boundary that is not 7. The symmetry properties of the cell problems 
are summarized in Table 5. We specify these Dirichlet boundary conditions on the veloc- 
ity together with the original boundary condition on 7. When forming the weak form of 
the problem, we use neutral boundary conditions, a • n = on the part of the boundary 
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that is not 7. For example, in solving the permeability cell problem of Section 3.1 for k 1 
and q±, the weak form is 

/ Vk 1 : V0 - giV • - • k 1 = / • d (Bl) 

351 where and ip are test functions. The Dirichlet boundary conditions are indicated in 

352 Figure 17. On the part of the boundary that is the interface, 7, we have applied the no - 
3b3 slip condition. 

B2. Solver Algorithms 

354 We discretize the Stokes equations for the cell problems using the P2-P1 formulation 

355 described in Elman et al. [2005]. The FEniCS libraries are used to generate code for the 

356 weak forms of the equations and assemble the associated matrices and vectors, [Dupont 

357 et al, 2003; Kirby and Logg, 2006, 2007; Logg, 2007]. These vectors and matrices are 

358 passed to PETSc and solved using algebraic Multigrid preconditioned GMRES, [Balay 

359 et al, 1997, 2004, 2001]. Domains and meshes were created with CUBIT [Sandia Corpo- 

360 ration, 2008]. The versions of the software we used are summarized in Table 6. 

To study problems with 0(10,000 - 100,000) elements and 0(100,000 - 1,000,000) 
unknowns, we rely on a Stokes preconditioner employing the pressure mass matrix of 
Elman et al. [2005]. The Stokes system is 







u 


= K 


u 




f" 


(b 


7) 


_p_ 


P_ 




s. 



where A is matrix corresponding to the weak form of the vector Laplacian, B T is the 
matrix corresponding to the weak form of the gradient, and B is the matrix corresponding 
to the weak form of the divergence. This is preconditioned with an approximate inverse 
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Of 




361 Q is the pressure-mass matrix. 

362 As our meshes are unstructured, the HYPRE library is used for algebraic multigrid 

363 preconditioning. In particular, we use BoomerAMG. We apply this on all of P, although 

364 we could have only used this on the A block, and relied on Jacobi or another light weight 

365 pre-conditioner for the Q block. 

B3. Examples and Benchmarks 

366 As a test, we solve the permeability cell problem of Section 3.1, with the symmetry 

367 reductions, for flow past a sphere of radius 0.3. It is meshed with a characteristic size of 
363 .03125, consisting of 119317 tetrahedrons. The results are summarized in Table 7 

369 For comparison, Jung and Torquato [2005] ran a time dependent problem to steady 

370 state and used an immersed boundary method with finite volumes. In our COMSOL 

371 computation, we used "fine" meshing, with 29649 elements, 134260 degrees of freedom, 

372 and a relative tolerance of le-10 in the solver. 

373 The objective function, (kj), converges as we refine our mesh; see Table 8 for a com- 

374 parison of different meshes for this problem. 

375 Table 9 summarizes the convergence results for flow around a sphere of radius 0.45. 

376 Again, our method appears to be quite effective. 

377 As another example, we solve the dilation stress cell problem from Section 3.2 . Solved 

378 on the domain complementary to a spherical inclusion of radius 0.2, the convergence 

379 results are summarized in Table 10. 
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380 The data in Tables 8- 10 were computed with default PETSc KSP tolerances. The 

381 automatically generated mesh was constructed with the CUBIT command: 

382 volume 6 sizing function type skeleton scale 3 time_accuracy_level 2 

383 min_size auto max_size 0.2 max_gradient 1.3 

384 While these convergence results are encouraging, our data is imperfect. Continuing 

385 with the dilation stress example, consider the data in Figure 18. Comparing Figures (a), 
sse (c), and (e), it would appear that the domains with smaller fluid inclusions have less 

387 well resolved pressure fields. They could likely be resolved with additional resolution. 

388 However, we use this data and believe it to be valid for several reasons: 

sag 1. It is preferable to have all domains meshed with the same algorithm. 

390 2. While the pressure fields may not be resolved, the error appears at the interface, 

391 and we are interested in the cell average. Moreover, the relative variations about the cell 

392 average are small. 

393 3. The corresponding velocity fields, with magnitudes pictured in Figures (b), (d), and 

394 (f), appear to be smooth, suggesting we are converging towards the analytical solution. 

395 4. The cell averages are consistent with the trends from the better resolved cases. 

B4. Cell Problem Data 

396 All meshes were generated using CUBIT with the command: 

397 volume 6 sizing function type skeleton scale 3 time_accuracy_level 2 

398 min_size auto max_size 0.2 max_gradient 1.3 

399 Problems were solved in PETSc with a relative tolerance of 1CT 8 and and absolute 

400 tolerance of 10 50 . 
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Table 1. Notation for macroscopic equations derived by homogenization. 

Symbol Meaning 

<5 C omp. Compaction length 

e(v) Strain rate tensor, e(v) = |(Vv + (Vv) T ) 

r/eff. Supplementary anisotropic viscosity derived by homogenization 

4> Porosity 

K Permeability tensor of the matrix derived by homogenization 

k e f{. Isotropic permeability of the matrix derived by homogenization 

jif Shear viscosity of the melt 

p s Shear viscosity of the matrix 

P Macroscopic (fluid) pressure derived by homogenization 

Pf Melt density 

p s Matrix density 

p Mean density, p = pf<p + (1 — (f>)p s 

V-f Macroscopic fluid velocity derived by homogenization 

V s Macroscopic solid velocity derived by homogenization 

Cefr. Bulk viscosity of the matrix derived by homogenization 
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Table 2. Notation for constitutive relations in other models. 



Symbol Meaning 

0* Critical porosity for activation of ji s+ f viscosity 

k Permeability of the matrix 

fi s+ f Shear viscosity of the matrix in the presence of melt 

( s Bulk viscosity of the matrix 



Table 3. Notation for cell problems. 

Symbol Meaning 

e Ratio of microscopic and macroscopic length scales, e = £/L 

e y (v) Strain rate tensor, e y {y) = |(V y v + (V y v) T ) 

7 Interface between melt and matrix within cell Y 

£ Grain length scale 

L Macroscopic length scale 

V y Gradient taken with respect to the y argument 

V y - Divergence taken with respect to the y argument 

Q Macroscopic region containing both melt and matrix 

y Coordinate within the cell, Y 

Y The unit cell 

Yf Portion of unit cell occupied by melt 

Y s Portion of unit cell occupied by matrix 

C Pressure of the cell problem for a unit forcing on the divergence equation 
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Figure 1. The macroscopic domain Q. The gray body is occupied by the matrix and 
the white inclusions are the melt. 




Figure 2. The cell domain, Y, divided into fluid and solid regions, Yf and Y s . The two 

phases meet on interface 7. 
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Table 4. Notation for cell problems, continued. 



Symbol Meaning 

X lm Velocity of the cell problem for a unit shear stress forcing on the solid in the 
Im component of the stress tensor 

ej Unit vector in the i-th coordinate, ef = (1, 0, 0) 

e ff. Effective porosity, portion of the porosity in which there is appreciable flow. 

4>cS. < <t>- 

k l Velocity of the cell problem for a unit forcing on the fluid in the ej direction 

\t\ First component of the velocity from the cell problem with unit forcing on the 
fluid in the ei direction 

£ Velocity of the cell problem for a unit forcing on the divergence equation 

(■) / Volume average of a quantity over the melt portion of a cell, (■) / = j Y -dy 

(■) s Volume average of a quantity over the matrix portion of a cell, (-) 8 = j Y -dy 

7r Pressure of the cell problem for a unit shear stress forcing on the solid in the 
lm component of the the stress tensor 

q l Pressure of the cell problem for a unit forcing on the fluid in the e« direction 

( Pressure of the cell problem for a unit forcing on the divergence equation 



Figure 3. A cell geometry composed of triply intersecting cylinders of equal radius. 




radius b 
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Porosity, <p 



Figure 4. Numerically computed values of (kj)j, the effective permeability, plotted 
against porosity for both the tube geometry and the sphere+tube geometry. The scattered 
circles are data from sphere+tube geometries, colored by the ratio of tube radius to sphere 
radius. The tube geometry offers an upper bound for a given porosity. 
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1(T 3 1(T 2 10- 

Porosity, </> 



Figure 5. Numerically computed data for the dilation stress cell problem, (9a-9c), on 
the three geometries, along with the least square fits (12), (26a), and (26b)). The scattered 
circles are data from sphere+tube geometries, colored by the ratio of tube radius to sphere 
radius. The sphere+tube data is bounded between the the tube data and the sphere data. 
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Porosity, 4> 



Figure 6. Numerically computed data for the normal stress cell problem on the three 
geometries, along with the least square fits (15), (27a), and (27b). The scattered circles 
are data from sphere+tube geometries, colored by the ratio of tube radius to sphere radius. 
At small porosity there is little variation amongst the simulated domains. 
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Figure 7. Numerically computed data for the shear stress cell problem on the three 
geometries, along with the least square fits (16), (28a), (28b). The scattered circles are 
data from sphere+tube geometries, colored by the ratio of tube radius to sphere radius. 
Sphere+Tube data points are constrained between the sphere data and the tube data. At 
all simulated porosities, there is less than an order of magnitude of variation. 
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melt channel 




Figure 8. An idealized Olivine grain from Figure 1 of Zhu and Hirth [2003]. Melt 
channels are found at triple junctions, while melt pockets are found quadruple junctions. 



radius b 




radius a 

Yf Y s 

Figure 9. A cell geometry composed of triply intersecting cylinders of equal radius, 
with a sphere at the intersection. 
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Figure 10. The magnitude of the velocity for a permeability cell problem. This 
corresponds to the sphere+tube domain with sphere radius a = .06 and tube radius 
b = .03. 




Figure 11. The magnitude of the velocity for a permeability cell problem. This 
corresponds to the sphere+tube domain with sphere radius a = .20 and tube radius 
b = .03. 
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Tube Radius 

Figure 12. (kj)/, the effective permeability, plotted against tube radius for both 
the tube geometry and the sphere+tube geometry. The scattered circles are data from 
sphere+tube geometries, colored by the ratio of tube radius to sphere radius. The tube 
geometry offers a lower bound for a given porosity. 
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Figure 13. The dark portion is the postulated effective porosity for the sphere+tube 
domains. 
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Effective Porosity, 4> e ^ 

Figure 14. (k})/, the permeability, plotted against the effective porosity for the 
sphere+tube geometry. The scattered circles are data from sphere+tube geometries, col- 
ored by the ratio of tube radius to sphere radius. The tube data plotted against porosity 
also appears. 
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(Effective Porosity) 1 88 (Porosity) 



Figure 15. (kj)/, the permeability, plotted against the effective porosity and porosity, 
using (23). The scattered circles are data from sphere+tube geometries, colored by the 
ratio of tube radius to sphere radius. The tube data is also plotted, substituting for 
4> eS . in (23). 
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Figure 16. Cell averaged pressure, (A7), plotted as a function of porosity. The data 
from the numerical solution of the cell problem given by equations (9a - 9c) and posed 
on y s cube also appears. It is in good agreement with the analytic solution on the spherical 
domain, Y s sphere . 
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Table 5. Summary of symmetry conditions used as Dirichlet boundary conditions in 

the cell problem computations. 



Cell Problem Velocity 


2/i = 0,- 


5 


2/2 = 


0,-.5 


2/3 = 0, 


-.5 


k 1 


^2 = ^3 = 





k\ 
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H = 







x\ l = o 




xf 


= 


xf = 





x 12 


X2 '1 X3 2 = 





x\ 2 = 


xl 2 = o 


xf = 







6 = 




6 
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6 = 
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symmetry reduced domain. 
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Table 6. Software versions 



Package 


Version 


CUBIT 


11.0 


DOLFIN(FEniCS) 


0.7.2 


FFC(FEniCS) 


0.4.4 


FIAT(FEniCS) 


0.3.4 


HYP RE 


2.0.0 


PETSc 


2.3.3 


UFC(FEniCS) 


1.1 


UMFPACK 


4.3 



Table 7. Convergence comparison between solvers 



Method 


(ki> 


BoomerAMG on P + GMRES 


0.0447051 


JT05 


0.045803 


COMSOL 


0.044497 
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Table 8. Convergence data for permeability cell problem I 



Mesh Size 


No. Cells 


No. d.o.f. 


(ki)/ 


|%A| 


0.25 


56 


430 


0.0492875 


— 


0.125 


534 


2950 


0.0472207 


0.0419336 


0.0625 


3673 


17988 


0.0452142 


0.042492 


0.03125 


26147 


119317 


0.0447051 


0.0112597 


auto 


14803 


70525 


0.0445419 





Table 9. Convergence data for permeability cell problem II 



Mesh Size 


No. Cells 


No. d.o.f. 


(ki)/ 


|%A| 


0.25 


61 


475 


0.00763404 




0.125 


384 


2302 


0.00651896 


0.146067 


0.0625 


2620 


13370 


0.00626809 


0.0384831 


0.03125 


19916 


93011 


0.00617889 


0.0142308 


auto 


23776 


112620 


0.00616139 




JT05 






0.0064803 




COMSOL 






0.006153 
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Table 10. Convergence data for dilation stress cell problem 



Mesh Size 


No. Cells 


No. d.o.f. 


(0, 


|%A| 


0.25 


73 


541 


46.6432 


— 


0.125 


514 


2862 


44.1747 


0.052923 


0.0625 


3813 


18575 


40.7177 


0.0782575 


0.03125 


29063 


132115 


39.4805 


0.0303848 


auto 


13725 


64205 


39.1558 




COMSOL 






39.117074 
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(e) (f) 



Figure 18. Figures on the left are the pressure fields for domains complementing 
spheres of radii a = .40, a = .20, and a = .10. Figures on the right are the corresponding 
velocity magnitude fields. 
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